home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / PROGRAMR / GSRC208A.ZIP / METCONAN.C < prev    next >
C/C++ Source or Header  |  1992-12-05  |  5KB  |  174 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*     metabolic control analysis    */
  11. /*        using Reder's method       */
  12. /*                                   */
  13. /*          MICROSOFT C 6.00         */
  14. /*           QuickC/WIN 1.0          */
  15. /*             ULTRIX cc             */
  16. /*              GNU gcc              */
  17. /*                                   */
  18. /*   (include here compilers that    */
  19. /*   compiled GEPASI successfully)   */
  20. /*                                   */
  21. /*************************************/
  22.  
  23.  
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <float.h>
  27. #include "globals.h"                       /* global symbols                 */
  28. #include "globvar.h"                       /* extern global variable         */
  29. #include "matrix.h"                        /* symbols & function declarations*/
  30. #include "datab.h"
  31.  
  32. #define MCA_OK 0
  33. #define MCA_NOCC 1
  34.  
  35. void calc_Dxv( void )
  36. {
  37.  register int i, j;
  38.  
  39.  for( i=0; i<nsteps; i++ )
  40.   for( j=0; j<nmetab; j++ )
  41.    Dxv[i][j] = partder[i][j]( xss, i, j );
  42. }
  43.  
  44. int calc_Gamma( void )
  45. {
  46.  register int k, j;
  47.  int i;
  48.  
  49.  for(i=0;i<indmet;i++)                               /* aux1 = rstoi * Dxv   */
  50.   for(j=0;j<nmetab;j++)
  51.   {
  52.    aux1[i][j] = (double) 0;
  53.    for(k=0;k<nsteps;k++) aux1[i][j] += (double) rstoi[i][k] * Dxv[k][j];
  54.   }
  55.  for(i=0;i<indmet;i++)                               /* aux2 = aux1 * ml     */
  56.   for(j=0;j<indmet;j++)
  57.   {
  58.    aux2[i][j] = (double) 0;
  59.    for(k=0;k<nmetab;k++) aux2[i][j] += aux1[i][k] * (double) ml[k][j];
  60.   }
  61.  if( lu_inverse( (double(*)[MAX_MET][MAX_MET])aux2,
  62.                  (double(*)[MAX_MET][MAX_MET])aux1,
  63.                  indmet
  64.                ) == NON_INVERT )                     /* aux1 = 1 / aux2      */
  65.  {
  66.   printf("  * non-invertible Dxv\n");
  67.   return MCA_NOCC;
  68.  }
  69.  for(i=0;i<nmetab;i++)                                 /* aux2 = ml * aux1     */
  70.   for(j=0;j<indmet;j++)
  71.   {
  72.    aux2[i][j] = (double) 0;
  73.    for(k=0;k<indmet;k++) aux2[i][j] += (double) ml[i][k] * aux1[k][j];
  74.   }
  75.  for(i=0;i<nmetab;i++)                               /* Gamma = -aux2*rstoi  */
  76.   for(j=0;j<nsteps;j++)
  77.   {
  78.    Gamma[i][j] = (double) 0;
  79.    for(k=0;k<indmet;k++) Gamma[i][j] -= aux2[i][k] * (double) rstoi[k][j];
  80.   }
  81.  return MCA_OK;
  82. }
  83.  
  84. void calc_C( void )
  85. {
  86.  register int i, j;
  87.  int k;
  88.  
  89.  for(i=0;i<nsteps;i++)                               /* aux = Dxv * Gamma    */
  90.   for(j=0;j<nsteps;j++)
  91.   {
  92.    aux1[i][j] = (double) 0;
  93.    for(k=0;k<nmetab;k++) aux1[i][j] += Dxv[i][k] * Gamma[k][j];
  94.   }
  95.  for(i=0;i<nsteps;i++)                               /* C = I + aux          */
  96.   for(j=0;j<nsteps;j++)
  97.    C[i][j] = aux1[i][j] + ( i==j ? (double) 1 : (double) 0 );
  98. }
  99.  
  100. void calc_tt( void )
  101. {
  102.  register int i,q;
  103.  
  104.  ttt = (double) 0;                                   /* reset total t.time   */
  105.  for( i=0; i<nmetab; i++)
  106.  {
  107.   tt[i] = (double) 0;                                /* reset trans. time    */
  108.   for( q=0; q<nsteps; q++ )
  109.    if ( stoi[i][q] > 0 )                             /* sum of positive fl.  */
  110.     tt[i] += fabs( rateq[q](xss, q) );
  111.   if ( tt[i] == (double) 0 )                         /* if no positive fl... */
  112.   {
  113.    for( q=0; q<nsteps; q++ )
  114.     if ( stoi[i][q] < 0 )                            /* sum of negative fl.  */
  115.      tt[i] += fabs( rateq[q](xss, q) );
  116.   }
  117.   if ( tt[i] != (double) 0 )
  118.   {
  119.    tt[i] = xss[i] / tt[i];
  120.    ttt += tt[i];
  121.   }
  122.   else
  123.   {
  124.    if (xss[i] == x0[i]) tt[i] = (double) 0;      /* disguised extern metablt */
  125.    if (xss[i] != x0[i]) ttt = tt[i] = INFINITY;  /* if rate=0 then tt=+oo    */
  126.   }
  127.  }
  128.  for( i=nmetab; i<totmet; i++ ) tt[i] = (double) 0;
  129. }
  130.  
  131. int equilibrium( void )
  132. {
  133.  register int i;
  134.  double hr;
  135.  
  136.  for( i=0, hr = (double) 0; i<nsteps; i++ )
  137.   if ( rateq[i](xss, i) > hr )
  138.    hr = rateq[i](xss, i);
  139.  
  140.  return ( hr <= options.hrcz );
  141. }
  142.  
  143. void scale_coef( void )
  144. {
  145.  int i, j;
  146.  
  147.  for( i=0; i<nsteps; i++ )
  148.   for(j=0;j<totmet;j++)
  149.   {
  150.    if( fabs( flux[i] ) >= options.hrcz )
  151.     Dxv[i][j] *= xss[j] / flux[i];
  152.    else Dxv[i][j] = FLT_MAX;
  153.   }
  154.  
  155.  for( i=0; i<totmet; i++)
  156.   for( j=0; j<nsteps; j++)
  157.   {
  158.    if( fabs( xss[i] ) >= options.hrcz )
  159.     Gamma[i][j] *= flux[j] / xss[i];
  160.    else
  161.     Gamma[i][j] = FLT_MAX;
  162.   }
  163.  
  164.  for( i=0; i<nsteps; i++)
  165.   for( j=0; j<nsteps; j++)
  166.   {
  167.    if(
  168.       ( fabs( flux[i] ) >= options.hrcz ) &&
  169.       ( fabs( flux[j] ) >= options.hrcz )
  170.      )
  171.     C[i][j] *= flux[j] / flux[i];
  172.    else C[i][j] = FLT_MAX;
  173.   }
  174. }